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Abstract 

We present extensive Monte-Carlo spin dynamics simulations of the clas- 
sical XY model in three dimensions on a simple cubic lattice with periodic 
boundary conditions. A recently developed efficient integration algorithm for 
the equations of motion is used, which allows a substantial improvement of 
statistics and large integration times. We find spin wave peaks in a wide 
range around the critical point and spin diffusion for all temperatures. At 
the critical point we find evidence for a violation of dynamic scaling in the 
sense that independent components of the dynamic structure factor S(q,u) 
require different dynamic exponents in order to obtain scaling. Below the 
critical point we investigate the dispersion relation of the spin waves and the 
linewidths of <S(q, u>) and find agreement with mode coupling theory. Apart 
from strong spin wave peaks we observe additional peaks in S*(q, uS) which can 
be attributed to two-spin wave interactions. The overall lineshapes are also 
discussed and compared to mode coupling predictions. Finally, we present first 
results for the transport coefficient D(q,u) of the out-of-plane magnetization 
component at the critical point, which is related to the thermal conductivity 
of 4 He near the superfluid-normal transition. 
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I. INTRODUCTION 



The theoretical investigation of classical spin systems has played a key role in the under- 
standing of phase transitions, critical behavior, scaling, and universality In particular, 
the claasical Ising, the XY, and the Heisenberg model are the most relevant spin models in 
three dimensions. Each of these simple models represents a universality class which, apart 
from the spatial dimensionality and the range of the interactions, is characterized by the 
number of components N of the order parameter, e.g, the magnetization in the case of 
ferromagnetic models. Despite their simplicity these spin systems continue to be of high 
relevance within the framework of dynamic behavior near critical points || (see also Ref. 

for a recent review). The Ising (N = 1), the XY (N = 2), and the Heisenberg (N = 3) 
universality class can be extended towards dynamic universality classes, which in addition 
to their static properties are characterized by the set of conservation laws 0. Special at- 
tention must be paid to the presence of energy conserving driving terms in the equations of 
motion which lead to propagating modes (spin waves) below the critical temperature and 
thus modify the dynamics || . The discrete nature of Ising spins does not allow such terms 
so that its dynamics is always of relaxational type (see Ref. || for a complete classification). 

The simplest spin model which allows propagating modes is a particular version of the 
ferromagnetic XY model (N = 2). The dynamics here is characterized by a nonconserved 
order parameter which is dynamically coupled to a conserved quantity (see Sec. II). The 
presence of spin waves reduces the value of the dynamic critical exponent z as compared to 
pure critical relaxation ||. The same is true for the isotropic Heisenberg model for which 
the N = 3 component magnetization vector is always conserved in the presence of energy 
conserving driving terms. If the model is ferromagnetic, the magnetization is the order pa- 
rameter. However, for an antiferromagnet, the magnetization plays the role of a conserved 
vector which is dynamically coupled to the nonconserved order parameter (staggered mag- 
netization). This difference in the conservation laws causes the classical Heisenberg ferro- 
and antiferromagnet to be in different dynamic universality classes although they belong to 
the same static universality class. Due to their fundamental role in the understanding of the 
critical dynamics in magnets Heisenberg ferro- and antiferromagnets have been thoroughly 
studied analytically by mode coupling theories (see Ref. || for a general overview) especially 
in presence of dipolar interactions |J and numerically by spin dynamics in d — 2 [J7| and in 
d = 3 HI! and by methods closely related to molecular dynamics jKj . 

The XY model may be viewed as a Heisenberg ferromagnet with an easy-plane (xy) 
anisotropy such that the order parameter has only two components. Planar ferromagnets 



are realized by layered compounds such as K 2 CuF 4 |nj and Rb 2 CrCl 4 |L2] which almost act 
as two-dimensional systems. The best results available today have been obtained on CoCl 2 
intercalated in graphite where a crossover from two-dimensional to three-dimensional 
behavior in the correlations has been observed below the Kosterlitz-Thouless temperature. 
Apart from the evident interpretation as a planar ferromagnet the XY model captures a 
larger variety of phenomena than the Ising or the Heisenberg model. Despite its continuous 
0(2) symmetry the XY model undergoes a continuous phase transition at a finite tempera- 



ture in two dimensions, known as the Kosterlitz - Thouless transition [14}|. Rather than by 



the onset of long-ranged order the transition is solely characterized by a diverging correlation 
length, when the critical temperature is approached from above. Due to a peculiar conspir- 
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acy between the spatial dimensionality d = 2 and the number of spin components N = 2 
configurations of bound and free vortices dominate the critical behavior of the XY model, 
where the unbinding of vortex pairs marks the point of the phase transition ||14|| . Naturally, 
many attempts to describe the critical dynamics of the XY model in d = 2 theoretically 
are based on the dynamics of vortices and vortex pairs fl5|-|T8|1 . According to analytical and 
numerical investigations for the ferromagnetic case [|16H T9|1 the in-plane component S xx (c[, uj) 
and the out-of-plane component S zz (q, uj) of the dynamic structure factor are expected to 
have central peaks above the transition. The line shapes of these peaks are predicted to be 
squared Lorentzian and Gaussian, respectively. Below the transition only spin-wave peaks 
are expected [16| . To test these specific predictions much numerical effort has been spent on 
spin dynamics simulations of the XY model in d = 2 Although dynamic finite-size 

scaling and the value of the dynamical exponent z = 1 has been confirmed to a high degree of 
confidence [pOR , the measured lineshapes of 5(q, uj) pOj^lj are apparently not well captured 



by analytical theories (see Ref. Q for details). It is therefore instructive to measure S^q, uj) 
for the XY model in d = 3 for which configurations of bound or free vortices do not play any 
special role for the critical behavior and should therefore not provide particularly noticable 
contributions to the structure factor. In d — 3 the dynamics of the planar ferromagnet 
has been investigated by mode coupling theory and specific predictions concerning line 
shapes and -widths have been made which can be compared with our data (see Sec. IV). 

It is well known that the A-transition of 4 He is the in the XY universality class, but 
the applications of the XY model for the physics of 4 He reach far beyond that. The spin 
dynamics for the XY model is the lattice analogue of the dynamical model E (symmetric 
planar ferromagnet M) which asymptotically also describes the critical dynamics of 4 He near 
the A-line |3,23-25|. If one therefore studies the transport properties of the XY model near 
the critical point T c , one should obtain lattice analogues of the corresponding transport co- 
efficients of 4 He near the A-transition. In this respect the aforementioned conserved quantity 
plays a particularly interesting part, because it is related to the entropy density in 4 He and 
its associated transport coefficient corresponds to the thermal conductivity of 4 He p3|j25 



which is an experimentally accessible quantity p6| . Below the critical temperature spin 
waves in the XY model then correspond to travelling waves of second sound in 4 He. These 
propagating modes cause the thermal conductivity to diverge at the lambda transition of 
bulk 4 He p3| , |25[| . In a finite system like our simulation sample, one therefore expects critical 
finite-size rounding of the thermal conductivity, which can be studied in the framework of 
the spin dynamics simulation and which should also be observable in experiments. 

To what extent the spin dynamics simulation actually captures the critical dynamics 
of 4 He is a rather delicate question. Although the asymptotic behavior is described by 
model E, the actual crossover to the asymptotic behavior, i.e., the decay of nonasymptotic 
corrections is governed by the specific heat exponent a ~ —0.013 f27| , which is so small that 
the true asymptotic behavior will never be seen in a simulation. From the point of view 
of analytic theory this means that in order to capture the nonansymptotic effects present 
in experiments one has to replace model E by the more complicated model F ^,23,2^ 



From the point of view of spin dynamics simulations this means that one has to look for 
sources of such nonasymptotic behavior artificially generated by the simulation method and 
other nonasymptotic corrections not captured by the model or the method (see Sec. II). 
Apart from these problems, it should also be mentioned that the dynamical model E has 
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two renormalization group fixed points. One of these fixed points yields dynamic scaling 
with a single dynamic exponent z = d/2 in d dimensions, whereas the other gives rise 
to a weak violation of dynamic scaling. Theoretical arguments f24] , |25| , |28|| and experimental 
evidence indicate that the latter fixed point is the stable one for 4 He in d = 3, i.e., the 
critical dynamics is characterized by two different dynamic exponents (order parameter) 
and z m (conserved quantity) which fulfill the scaling relation z^ + z m = d. Their difference 
u w = — z m ^ has the nature of a dynamic Wegner exponent and is known as the 



transient exponent 25 - 2 



The remainder of the paper is organized as follows. In Sec. II we present the model and 
the simulation methods used to generate equilibrium configurations and to obtain the critical 
point of the model and its static critical exponents. Furthermore, the equations of motion 
and the method used to integrate them numerically are presented. In Sec. Ill we briefly 
discuss the static critical behavior of our model and present an accurate estimate of the 
critical temperature. Sec. IV is devoted to the discussion of the dynamic structure factor and 
the comparison with predictions of analytic theory. In Sec.V we present first results for the 
lattice analogue of the thermal conductivity and discuss its scaling properties. A summary 
and prospects for future work are given in Sec. VI. Unless otherwise stated statistical errors 
quoted in this work correspond to one standard deviation. 



II. MODEL AND SIMULATION METHOD 

The system under investigation is given by a ferromagnetic Heisenberg model with the 
strongest possible easy-plane anisotropy. The model Hamiltonian reads 

H = -JJ2 ( S ? S j + S i S j) > (2- 1 ) 

(ij) 

where (ij) denotes a nearest neighbor pair of spins on a simple cubic lattice in three dimen- 
sions. The lattice contains L lattice sites in each direction and in order to avoid surface effects 
periodic boundary conditions are applied. Each spin Sj is a classical spin Sj = (S?, Sf , S-) 
with the normalization |Sj| = 1. The easy-plane anisotropy in Eq.( |2.1| ) is the strongest pos- 
sible in the sense that the z-components of the spins do not couple, so that Eq. (|2.1| ) looks 
like the standard Hamiltonian for the usual (plane rotator) XY model. 

As a starting point for the spin dynamics a sequence of equilibrium configurations is 
needed to provide initial conditions for the equations of motion. These configurations are 
obtained from a Monte-Carlo simulation of the model Hamiltonian given by Eq.( |2.1| ). The 
Monte-Carlo algorithm chosen is a hybrid scheme, where each hybrid Monte-Carlo step 
(MCS) consists of 10 updates each of which can be one of: one Metropolis sweep of the 
whole lattice, one single cluster Wolff update f2~9fl , or one over relaxation update of the 



whole lattice || . The Metropolis algorithm updates the lattice sequentially in the standard 
way. According to the detailed balance condition we choose the acceptance probability 
p(f3AE) = 1 / (exp(j3 AE) + 1) for a single spin flip, where AE is the change in configurational 
energy according to Eq. (|2.1|) and (3 = 1/(/cbT). 



The Wolff algorithm also works the standard way [23], except that only the x and y 
components of the spins are used for the cluster growth. This means that a cluster update 
never changes the z-component of any spin so that the Wolff algorithm is nonergodic in 
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this case. Our cluster update is still a valid Monte-Carlo step in the sense that it fulfills 
detailed balance, however, in order to provide a valid Monte-Carlo algorithm it has to be 
used together with the Metropolis algorithm described above in a hybrid fashion. The 
use of Wolff updates allows us to take advantage of improved estimators for magnetic 
quantities. 

The overrelaxation part of the algorithm performs a microcanonical update of the con- 
figuration in the following way. The local configurational energy has the functional form of 
a scalar product of the spins, where according to Eq. (|2.1| ) only the x and y components are 
involved. With respect to the sum of its nearest neighbor spins each spin has a transverse 
component in the xy plane which does not enter the scalar product. The overrelaxation al- 
gorithm simply scans the lattice sequentially, determines this transverse component for each 
lattice site and flips its sign. This does not change the local configurational energy (AE = 0) 
and by virtue of the usual Metropolis acceptance function f(f3AE) = min(exp(— /3AE), 1) 
the update is always accepted. Along with this simple operation the sign of S- is flipped 
with probability 1/2 at each lattice site which according to Eq.( [2.1|) also does not change 
the energy of the configuration. This overrelaxation algorithm is similar to the one used 
in Ref. [§] and it quite efficiently decorrelates subsequent configurations over a wider range 
of temperatures around the critical point than does the Wolff algorithm. Typically, we use 
three Metropolis (M), five single cluster Wolff (C), and two overrelaxation updates (O) in a 
hybrid Monte-Carlo step in the critical region of our XY model. The inividual updates are 
mixed automatically in the program so that the update sequence (MCCMOCMCCO) 
is generated as one hybrid Monte-Carlo step in this case. The random number generator we 
use is the shift register generator R1279 given by the recursion relation X n = X„_ p © X n _ g 
for (p,q) = (1279, 1063). Generators like this are known to cause systematic errors in com- 
bination with the Wolff algorithm pi| ; however, for lags (p, q) as large as the ones used here 
these errors will be far smaller than typical statistical errors. They are further reduced by 
the hybrid nature of our algorithm | |32| . 

The spin dynamics of the XY model is defined by the equations of motion 



where 7i ist the Hamiltonian defined by Eq.( |2.1| ). One may interpret Eq. Q2.2p as the direct 



classical analogue of the Heisenberg equations of motion for spin operators, where h = 1 
so that energies and frequencies are measured in the same units. From the symmetry of 
Eq.( |2.1| ) it is evident that the components M x and M y of the magnetization M = J2k&k 
are not conserved under the dynamics given by Eq. (|2.2| ). Note that the two-component 
vector (M x , M y ) is the order parameter of the XY model. The z- or out-of-plane component 
M z of the magnetization M is just the conserved quantity within the framework of model 
E dynamics we have already referred to in Sec. I 0]. Note that Eq. Q2.ip is invariant with 
respect to the transformation M z — > —M z which is a symmetry required by model E. 

For the comparison of the critical spin dynamics with the critical dynamics according to 
model E it is important to realize, that the configurational energy is an additional constant 
of motion, because Eq. (|2~^ ), in contrast to the coarse grained model E, does not contain 
relaxation. Whether energy conservation is a reasonable assumption for the dynamics of the 
XY model or any other classical spin model is a question of the time scales to be resolved. 
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The most important time scale for our investigation is set by the propagating modes (spin 
waves) in the system and for these the configurational energy is indeed constant. Within 
the time scale of the spin waves thermal averages can therefore be replaced by averages over 
the initial configurations from which the time integration of Eq.( [2.2|) is started. For much 
longer times relaxation processes (equilibration with the heat bath) come into play which 
violate energy conservation and render our spin dynamics approach invalid. In the vicinity of 
the critical point energy conservation becomes particularly important, because the dynamic 
universality class may change under the influence of an additional conservation law. If model 
E is augmented by energy conservation (model E', see Ref. [Q) it turns out that the energy 
asymptotically decouples from the order parameter (M x , M y ) and the conserved out-of-plane 
magnetization M z and model E critical behavior is restored. However, energy conservation 
may introduce corrections to the asymptotic finite-size scaling behavior []23| which decay 
very slowly and may cause ambiguities in the scaling analysis of the spin dynamics data. 
Note that these corrections are generated by the spin dynamics method. 

The equations of motion given by Eq. (|2.2|) are integrated numerically for each initial spin 
configuration by a recently developed decomposition method |33] . This method guarantees 
exact energy conservation and conservation of spin length |Sfe| = 1 and conserves M z within 
its numerical truncation errors. For the present study a second order integrator is used with 
the time step 5t = 0.05/ J. This time step guarantees sufficient accuracy with respect to the 
conservation of M z and is much faster than well-known predictor corrector methods Pl,|33 



For some accuracy and stability tests the time step has been increased to 5t = 0.1/ J which 
still yields sufficient accuracy for the dynamic structure factor. Fourth order integrators are 
much more accurate as far as M z conservation is concerned, but their internal complexity 
makes them much slower than a second order method for the same time step | 33fl . Moreover, 
statistical errors are not decrased significantly by fourth order methods and we therefore only 
report results obtained by the second order method. The equations of motion are integrated 
to a final time of 800/ J and thermal averages are taken over 1000 initial configurations. All 
error bars of dynamic quantities correspond to one standard deviation. The simulations have 
been performed on various DEC alpha AXP, IBM RS6000, and HP RISC8000 workstations 
both at the RWTH Aachen and the BUGH Wuppertal. 



III. STATIC PROPERTIES OF THE XY MODEL 

A. Thermodynamic properties 

The basic ingredient for the spin dynamics simulation is provided by the sequence of 
initial spin configurations, which has to be generated according to the canonical ensemble 
in order to provide well defined thermal averages. Therefore, the static behavior of the XY 
model and especially the location of T c have to be determined first. For this purpose we 
employ the hybrid Monte-Carlo scheme described above for lattice sizes L between L = 20 
and L = 80. For each system size and temperature we perform 10 blocks of 10 3 hybrid steps 
for equilibration followed by 10 4 hybrid steps for measurements. Each measurement block 
yields an estimate for all static quantities of interest and from these we obtain our final esti- 
mates and estimates of their statistical error following standard procedures. The integrated 
autocorrelation time of our hybrid algorithm is determined by the autocorrelation function 
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of the energy or, equivalently, the modulus y + of the order parameter, which yield 
the slowest modes for the Wolff algorithm. The autocorrelation times are generally rather 
short, at T c (see below) they range from about 5 hybrid MCS for L = 20 to about 10 hybrid 
MCS for L = 80. For comparison, the autocorrelation function for the order parameter itself 
yields an autocorrelation time of less than one hybrid MCS. The values for the equilibration 
and measurement periods given above thus translate to roughly 100 and 1000 autocorrela- 
tion times, respectively, which is sufficient for all practical purposes. In order to obtain the 
best statistics for magnetic quantities a measurement is made after every hybrid MCS. 

The critical temperature T c is determined by monitoring the temperature and size de- 
pendence of the Binder cumulant ratio. Specifically, we measure the cumulants 

. (Mi) , , (jM + mi?) 

It turns out that u\ is more sensitive to changes in temperature and system size so that we 
only use u\ for the final fine-tuning of the temperature T. For convenience the temperature 
is expressed as the dimensionless reduced coupling K = J/(ksT), where K c = J/(/csT c ) 
denotes the critical point. From standard procedures []33| we obtain K c = 0.64440 ±0.00005. 



For comparison we mention that K c = 0.45420 ± 0.00002 for the standard plane rotator XY 
model on a simple cubic lattice in d = 3 [ . By ignoring corrections to scaling and averaging 



over the measurements for L = 40, 50, 60, and 80 we find the estimates 

u\ = 0.3789 ±0.0015 and u* 2 = 0.5859 ± 0.0008 (3.2) 

for the values u\ and u\ of the cumulants defined by Eq. fl3.lfl at the critical point. Note that 
for K = 0.6444 the value of u% remains within two standard deviations of u\ for all L. Our 
estimate for u 2 is within two standard deviations of the corresponding estimate 0.5891 ± 
0.0020 found in Ref. P5|] which already gives some evidence that the planar Heisenberg 



variant of the XY model studied here is indeed a member of the static XY universality class. 

The critical exponents are estimated from the critical finite-size scaling behavior of the 
average modulus {\JM% + M^j of the order parameter, the average square (Mj ± My) of 
the order parameter, and the temperature derivative of the latter. At T = T c , i.e., K = K, 
one finds the leading scaling behavior 



L-' 3 (Ml + Mf) ~ L 7/V , (3.3) 



y 

^\n{Ml + Ml) -I}* 

with the system size L, where (3, 7, and v are the critical exponents of the order parameter, 
the susceptibility, and the correlation length, respectively. During the data analysis it turns 
out that corrections to scaling can be ignored within the statistical error of the quantities 
in Eq.flpD. From our estimate K c = 0.6444 and L = 20, 24, 30, 36, 40, 50, 60, and 80 we 
find the following values for the critical exponents: 

P/v = 0.5179 ±0.0024, 

7/1/ = 1.965 ±0.005, (3.4) 
l/v = 1.494 ±0.013. 
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These values satisfy the scaling relation 2/3 /u + j/u = d. Furthermore, we directly obtain 
from Eq.(l3l) v = 0.6693 ± 0.0058, rj = 0.035 ± 0.005, 7 = 1.315 ± 0.012, and (3 = 0.3467 ± 
0.0034 in very good agreement with previous Monte-Carlo estimates [|35] and renormalization 
group theory for the 0(N = 2) Ginzburg-Landau model |36|. The critical exponent a of the 



specific heat can be obtained from, e.g. the hyperscaling relation, but the statistical error 
of v is too large to exclude logarithmic behavior (a = 0). Apart from this deficiency our 
simulations confirm XY-like critical behavior for our version of the XY model (see Eq. fl2.lp ) 
quite accurately. In the following k^T is measured in units of J chosen such that /csT c = 1. 



B. Static structure factor at criticality 

We continue this section with a short discussion of the static spin-spin correlation function 
(structure factor) G(q) at the critical temperature. The static structure factor is the spatial 
Fourier transform of the spin-spin correlation function 

Gaf> (Hi - Rj) = (S?Sl) - (S?)(Sf), (3.5) 

where a, (3 refer to spin the components, Rj and Rj denote lattice vectors, and the thermal 
average is indicated by (...). Note that Q a p = Q a p (Rj — Rj) = Q a p (Rj — Rj) due to trans- 
lational invariance and reflection symmetry of the system. The spatial Fourier transform is 
therefore given by a cosine transform which we use in the form 

GWK<l) = H Gap (Ri) cos q ■ Ri. (3.6) 

i 

With respect to the spontaneous magnetization of the XY model G a p(<\) has a transverse 
component Gj(q) and a longitudinal component Gj(q). In our Monte-Carlo simulation we 
measure G t and Gi by rotating the coordinate system of the spins around the z-axis such that 
the random but finite magnetization vector (M x , M y ) is aligned with the ^/-direction. The x- 
and ^-components of the spins in the rotated frame then correspond to the transverse and the 
longitudinal spin components, respectively, and their correlation functions yield G t = G xx 
and Gi = G yy . It should also be mentioned that the out-of-plane component G zz (q) of 
the static structure factor is independent of q and does not show any critical behavior. 
According to Eq. (|5TB|) the normalization of G a p(q) is such that G Q/ g(q = 0) = ksTxa/3, 
where Xa/3 is the static susceptibility. Note that G*(q = 0) = by definition and that 
Gi(q = 0) = ksTx' 1 where x' is the magnetic susceptibility with respect to the modulus of 
the magnetization. We furthermore limit the discussion to the (100) direction, because other 
lattice directions do not provide new information on a simple cubic lattice. The components 
of q are always measured in units of the inverse lattice contant. 

At T c the longitudinal component Gi(q) of the static structure factor can be described 
by the model function 

G l (q=(q,0,0)) = L 2 -Vg l (qL)h(q), (3.7) 
where 2 — r] = 7/1/ is taken from Eq.( J3.4j ), h(q) = [(g/2)/sin(g/2)] 2 captures lattice ef- 



fects [57], and the finite-size scaling function gi(x) is chosen as a simple generalization of a 



Lorentzian: 
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gi ( x ) = [l + (x/b) 2 } (2 "' V2 . (3.8) 

The parameters a and b are determined from a fit to the data, where a and b are independent 
of L. We have chosen a = 3.70 and b = 4.35 as obtained from a fit to the data for L = 60. 
For L = 40 and L = 80 a and 6 are found within less than 1% of these values, for L = 20 they 
are about 3% smaller. Despite its simplicity the model function given by Eq.fl3.7p captures 
the shape of the Monte-Carlo estimate for the longitudinal component of the structure factor 
remarkably well. The finite-size scaling analysis of our data for G/(q) for small q, is shown 
in Fig.|TJ. Deviations from finite-size scaling set in at qL = Air for L = 20, the data for 
L = 40, 60, and 80 collapse within the error bars up to qL = 6tt. Within the symbol sizes 
data collapse is obtained up to qL = 10ir, where lattice effects set in. For comparison the 
model scaling function gi(qL) (a = 3.70 and b = 4.35, see Eq. (|3.8fj ) is shown by the dashed 
line in Fig.|I|. The true scaling form of Gi(q) for small q is captured rather well by gi(x). 
However, the choice of the model function is, of course, not unique. 



IV. THE DYNAMIC STRUCTURE FACTOR 

The dynamic structure factor S(q, u) is the space-time Fourier transform of the position 
and time displaced spin-spin correlation function S a p (Rj — Rj, \t — t'\) which we define by 
(see Eq.Q) 

S aP (R, - R„ |t - 1'|) = (S?(t)Sf (t')> - (S?(t))(Sf(t')). (4.1) 

The indices a, (5 refer to spin components, Rj and R., are lattice vectors, and t and t' are mo- 
ments in time to which the initial spin configuration has evolved according to the equations 
of motion (see Eq. fl2.2| )). The average (...) is taken over the set of initial configurations as 
described in Sec. II. Note that S a p is also symmetric with respect to Rj and Rj (see Eq. ( |3.5|) ). 
As the second argument only the time displacement enters, because the equations of motion 
are invariant under the transformation Sj — ► —Si, t — > —t (see Eqs. (|2.1| ) and ( |2.2j )). The 
space-time Fourier transform is therefore also given by the cosine transform (see Eq. fl3.6Q ) 

POO 

Sosfq, a>) = 2/7T / S al 3 (Rj, Itl) cosq-Rj cosujt dt. (4.2) 
Jo - 

Note that the normalization in Eq.(|4.2|) is such that J °° 5 Q/ a(q, Lu)du = G aj g(q), where G Qj g(q) 
is the static structure factor defined by Eq. flSTB"!) . 

The out-of-plane component S zz (q,u) of the dynamic structure factor is associated with 
the conserved out-of-plane magnetization M z , i.e., ^ «S M (Rj, |t|) (see also Eq.(P~2|)) is a 
constant in time. The in-plane magnetization (order parameter) (M x ,M y ) is not conserved 
under the spin dynamics. Although the time scale set by the motion of M x and M y is 
considerably larger than typical time scales set by spin waves, all initial differences between 
longitudinal and transverse components of the spin correlations completely disappear during 
the integration of the equations of motion. We therefore only discuss the average of the xx 
and the yy component of S(q,u) and refer to it as the "in-plane component" S xx (q,uj). 

We now turn to the discussion of £>(q, to) above, below, and at the critical point. The 
correlation functions are measured up to a time displacement of 400/ J away from critical- 
ity, where only smaller systems are considered. At the critical temperature we measure 
correlations up to 600/ J for systems with up to L = 60 lattice sites in each direction. 
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A. The structure factor above T c 



LUQ, = ; , r/ Z ?Z , , l2 + - u T' tT , , ^ , (4-4) 



In order to avoid effects of criticality we choose the temperature T = 1.1T C in the 
following. Due to the absence of critical finite-size scaling at this temperature we limit 
the spin dynamics simulations to smaller systems with L = 20, 24, and 30. For better 
momentum resolution we only present results obtained for L = 30, smaller L yield identical 
results with a lower resolution in q. We also limit the presentation to S(q,u) in the (100) 
direction, other lattice directions provide essentially the same information. In Fig.0 we show 
S xx (q = (<?, 0, 0), oj) and S xx (q = (q, 0, 0),uj) for q = n/15 as functions of the dimensionless 
frequency u/ J. Both components apparently display a central peak without any additional 
features. For the values of T and q used in Fig.[| one expects a central peak of Lorentzian 
shape |f22| . As displayed by the solid line in Fig.^ a simple Lorentzian of the form 

L„(q,u) = fj^\ (4.3) 

characterized by an amplitude A xx (q) and a width T xx (q) captures the shape of S xx very 
well up to uj/J ~ 0.4, where the intensity has already dropped by an order of magnitude. 
For larger uj S xx decays faster than a Lorentzian. If one tries to fit the line shape of S zz also 
with a Lorentzian for small uj, it turns out that a better fit is obtained from a superposition 
of symmetrically placed Lorentzians. We use 

AM | A zz (q) 

l + [(u- uj(q))/T zz (q)} 2 1 + [(uj + uj(q))/T zz (q)f 

where uj{q) denotes the spin wave frequency which is used as an additional fit parameter. 
The dashed line in FigJ^ displays the Lorentzian fit for S zz according to Eq. (|4.4|) , where 
u;(7r/15) ~ 0.030 is indeed finite within its statistical error of roughly 10~ 3 . However, 
Eq.flPD only captures the shape of S zz up to u/J ~ 0.15. If Eq-flPQ (^(V 15 ) = °) is 
used instead, this frequency range becomes even smaller. This discrepancy in frequency 
range between the in-plane and the out-of-plane components of S"(q, uj) may be due to the 
difference in time scales between in-plane and out-of-plane modes as predicted by mode 
coupling theory p2 j. 

Spin wave signatures are still visible in S zz at T = 1.1T C as shown in Fig.^|, where S zz 
is plotted for the first four momenta q = mr/15, n = 1,2,3,4. For n > 2 a very broad 
maximum becomes visible which moves to the right as q increases. The appearance of spin 
wave signatures in S zz above T c is expected from renormalization group theory for model E 
dynamics P5] . For the q values used in Fig.|3] the shape of S zz is captured quite well by the 
Lorentzian defined by Eq.( j4.4|) . Due to the small enhancement of the line intensity over the 
signal level at uj = the peak position uj(q) and linewidth T zz (q) are not very well defined. 
We therefore refrain from a more detailed analysis at this point. The frequency dependence 
of the in-plane component S xx is dominated by a central peak as in FigfJ for all q. Only for 
large q near the Brillouin zone boundary a shoulder appears in S xx near the position of the 
spin- wave peak in S zz . We illustrate this in Fig.[|, where S xx and S zz are shown for q = n. 
Note that S zz , which is much smaller than S xx for small q (see Fig.[|), becomes comparable 
to S xx in magnitude for large q. In a qualitative sense we have recovered the same behavior 
as observed by spin dynamics simulations of the XY model in d = 2 above the Kosterlitz - 
Thouless transition (2(| . 
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B. The structure factor below T c 



As before we want to avoid critical finite-size effects also below T c and we therefore 
choose T = 0.9T C in the following. Again relatively small systems are sufficient for this 
investigation. Our main results have been obtained for L = 30, smaller systems yield the 
same results at a lower q resolution. The dynamics of the XY model below T c is dominated 
by spin waves which are visible in 5(q, u) as pronounced peaks at the spin wave frequency 
uj(q) as shown in Fig.|5| for q = 7r/15 in the (100) direction. The peak intensity of S xx is more 
than one order of magnitude larger than the corresponding peak intensity of S zz . S xx also 
displays a pronounced central peak, where the intensity exceeds the intensity of S zz by over 
three orders of magnitude. In a qualitative sense this is again the same behavior as observed 
in d — 2 ||20|| , where the central peak in S xx was not expected by analytical theory [TO. We 
also find some 'fine structure' at low intensities in S xx for which no theoretical predictions 
exist. In Fig.|5] we have marked these additional resonances by arrows. As demontrated 
already in d = 2 [^0J these signals can be interpreted as two-spin wave peaks in the following 



way. In order to produce a contribution to S xx at q = ir/15 in the (100) direction one 
can combine two spin waves at q = 7r/15, one in the (010) direction, which is eqivalent to 
the (100) direction on a simple cubic lattice, and one in the (110) direction (not shown), 
where the latter has a higher energy (frequency). The difference between the corresponding 
frequencies is marked as (1) in Fig.[5], their sum is marked as (2). A second way to get a 
contribution to S xx at q = 7r/15 in the (100) direction is to combine two spin waves in the 
(100) direction, one at q = it/ 15 the other at q = 27r/15, where again the latter has a higher 
energy. The sum of the corresponding frequencies is marked as (3) in Fig.[5] their difference 
almost coincides with the position of the spin wave peak and can therefore not be resolved. 

In order to monitor the q dependence of the spin wave frequency uj(q) and the line widths 
F xx (q) and T zz {q) of the spin wave peaks, we again employ fits to simple Lorentzians. For 
the in-plane component we generalize Eq . (|4 . 3|) to include the spin wave peaks: 



L ( a LU )= + A **^ + ^Wg) u 5) 

XXK ' 1 + l + [(cv-u;(q))/T xx (q)] 2 1 + [(« + u(q))/T xx (q)] 2 ' 



For the out-of-plane component we use Eq.( [4.4| ) which captures the shape of S zz for suffi- 
ciently small frequencies in a satisfactory way. The dispersion relation oj{q) can be obtained 
from Eq.( [4.5| ) or Eq. fl4.4j ), where the latter yields slightly smaller error bars, because S zz 
displays a sharper spin wave maximum. The estimates for u>(q), T xx (q), and T zz (q) depend 
on the frequency range over which Eqs.([4.5D and ( |4.4j ) are fitted to S xx and S zz , respectively. 
We have chosen a frequency window around the spin wave peak, where for S xx the central 
peak has been subtracted first. The error in the dispersion relation and the line widths 
is estimated by varying the size of the frequency window from about 1.2 times the half 
width to about twice the half width of the peak. The dispersion relation u)(q) obtained from 
this procedure for the (100) direction is shown in Fig.|6], the correponding zero temperature 
dispersion relation shown by the solid line. According to linear spin wave theory one obtains 

u(q)/J = 2V2dsin(q/2), (4.6) 

for T = where d — 3 in our case. As expected, the spin wave frequencies are 'renormalized' 
to fall below the T = dispersion curve. The functional form of ou(q) can be captured by a 
Fourier series. A convenient choice is 
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u (q) / J = a>\ sin(g/2) + a 2 sin(3g/2) + a 3 sin(5g/2) + . . . , 



(4.7) 



where the Fourier coefficients ai,a 2 ,a 3 , ■ ■ ■ are rapidly decreasing. If Eq. (|4.7j ) is truncated 
after the second term one obtains a x = 2.919 and a 2 = —0.116 from a least square fit which 
connects all data points within their error bars as shown by the dashed line in Fig|| 

The line width T xx (q) of the spin wave peak in the in-plane component of the structure 
factor is shown in Fig.|7|(a). As before the data can be analyzed by a Fourier series with 
rapidly decreasing coefficients. A convenient choice here is 



where a refers to x (in-plane) and z (out-of-plane). A fit to T xx with only two coefficients 
yields bo = 0.4622 and b\ = —0.4559 which is shown by the dashed line in Fig.^(a). The 
statistical error of bo and &i is about 5 x 10~ 4 . The data for small q are represented quite 
well, whereas near the Brillouin zone boundary significant deviations occur. These can be 
reduced very quickly by including higher Fourier modes in the fit (not shown), where the 
higher Fourier coefficients decrease rapidly in magnitude. Here, we are primarily interested 
in the q dependence of the line width near the center of the Brillouin zone (see below). 
A corresponding analysis has been performed for the line width T zz (q) of the spin wave 
peak in the out-of-plane component of the structure factor. The result is shown in Fig. 0(b), 
where the Fourier coefficients in the fit (dashed line) are given by bo = 0.4160(5) and 
b\ = —0.4165(5). Again the small q behavior is captured very well by the fit, but near 
the Brillouin zone boundary deviations occur which can also be reduced very quickly by 
including higher Fourier modes. 

The limit q — > is of particular interest for the line widths, because it reflects the 
influence of conservation laws on the dynamics. The line width T xx (q = 0) can be interpreted 
as the relaxation rate of the (non-conserved) order parameter and therefore T xx (q = 0) 
should be positive. The line width T zz (q = 0) is the relaxation rate of the out-of-plane 
magnetization M z which is conserved, i.e., T zz (q = 0) should vanish. From the Fourier fits 
shown in Fig.0 one finds T xx (0) ~ 0.0063 and r„(0) ~ -0.0005 with statistical errors of 
about 7 x 10~ 4 in both cases. The extrapolation of T xx (q) and T zz (q) to q = is shown in 
Fig.||]. The q dependence of the line widths for small q is quadratic as anticipated by the 
Fourier fits shown in FigfJ. From a fit to a straight line we find T xx (0) = 0.0076±0.0006 > 
(solid line), whereas T zz (0) = (dashed line) within its statistical error in agreement with 
the conservation laws. Finally, we note that the q 2 dependence of the linewidth T zz is in 
agreement with the prediction of mode coupling theory |22| . 



The exploration of critical dynamics and dynamic scaling naturally requires most of the 
numerical effort. In order to reach the scaling regime large systems are required and we have 
therefore performed simulations for L = 20, 24, 30, 40, and 60. Our prime objective here is 
the test of dynamic finite-size scaling which we assume to be valid in the form |8|j20[| 



^aa(q)/J = b + h cosg + b 2 cos2g + . . . , 



(4.8) 



C. The structure factor at T c 



S aa (q,co)/G aa (q) = L Za E 



a a 



(qL,uL z «) 



(4.9) 
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at the critical point, where G aa is the static structure factor discussed in Sec. Ill and z a 
denotes the dynamic critical exponent. Note, that the q dependence is reduced to a depen- 
dence only ong= |q|, i.e., isotropic scaling in space has been assumed. The index a again 
refers to the spin component x (in-plane) or z (out-of-plane), where z x = z^ and z z = z m , 
respectively. In order to estimate z a we scale our data according to Eq.([4.9|) and test the 
result for data collapse in the limit of small frequencies. The static structure factor needed 
for normalization in Eq.( (4.9| ) is given by the equal-time correlations which we also extract 
from the spin dynamics data for consistency. 

In order to provide numerical estimates of the dynamic exponents z m and z^ first we 
analyze the dispersion relations w zz (q) and oo xx (q) at the critical point as obtained from a 
Lorentz fit to the spin wave peak of S zz (q,u) (see Eq. (|4.4j )) and a Lorentz fit to S xx (q, uj) 
according to Eq. ( |4.3| ) . The line shape of S xx is dominated by a strong central peak (see 
below) so that we restrict the analysis of the full dispersion relation to u zz (q). The result in 
the (100) direction is shown in Fig.|^ for L = 24 and L = 40. The data apparently collapse 
onto a single curve and show a linear behavior for q < n/2 down to q ~ v/20, where uj(q) 
becomes nonlinear. From finite-size scaling one expects the scaling form 

u aa = q Za Q a (qL) (4.10) 

at T c for sufficiently small q. For qL = 2ir and L = 20, 24, 30, 40 and 60 we obtain an 
approximation of the small q behavior of uo zz which is shown in Fig.[l(](a). According to 



Eq.( f4.10|) oo zz should vary as q Zm for fixed qL. A least square fit to the data shown in 



Fig.|T0"|(a) yields the dynamic exponent 

z m = 1.38 ±0.05 (4.11) 



which differs substantially from the mode-coupling prediction z m = z^ = 1.5 |[22| . From 
the scaling relation z^ = 3 — z m [p4] , p5| , p8| we obtain z^ = 1.62 ± 0.05 which can be tested 
against the small q behavior of u xx (q). The result is shown in Fig.|l0|(b). For the larger 
systems L = 30, 40, and 60 the data agree with the power law, but for smaller lattice 
sizes L = 20 and L = 24 systematic deviations occur. If we exclude these smaller systems 
from the power-law fit shown in Fig.|T0|(a) we find z m = 1.43 ± 0.14 which is at the upper 
error boundary of the previous estimate given by Eq.( 4.1l| ). An alternative estimate of the 



J ota 



dynamic exponents is provided by the median frequency tu^ a which is defined by the relation 

S aa (q, uo)duJ = |G aQ (q). (4.12) 

Note that we have normalized the dynamic structure factor such that its integral over all 
frequencies yields the static structure factor. We evaluate Eq.( |4.12j ) with the trapezoidal 
rule. Systematic errors induced by this simple numerical integration method are smaller 
than statistical uncertainties coming from the statistical errors of S xx (q,u) and S zz (q,uj). 
If the integration is performed with Simpsons rule the same results are obtained within 
statistical errors. According to Eq. (|4.9|) we expect the scaling behavior 

u aa = L~ Za M aa (qL) (4.13) 

of the median frequency at the critical point. Due to the presence of a strong central peak 
in S xx (q,u) the median frequency uj^ x for L = 60 turns out to be of the same size as the 
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frequency resolution of our data. We therefore limit the analysis of u™ x to the system sizes 
L = 20, 24, 30, and 40. The other median frequency io™ z can be determined accurately for 



all system sizes. The result for qL = 2tc is displayed in Fig.[Tl|. From a least square fit of 
u xx t° a power law for L = 20, 24, 30, and 40 we obtain z$ = 1.61 ± 0.03. If only the 
data for L = 30 and 40 are used we obtain = 1.65 ± 0.10. Both estimates are consistent 
with Eq. (|4.11|) and the scaling relation z ( p + z m = 3. The latter estimate is displayed by the 
solid line in Fig.|TT[ The corresponding fit of for L > 30 yields z m = 1.35 ± 0.01 which 
is also within the error bar of our previous estimate (see Eq.( [4.1lD , dashed line in Fig.|i~l"D. 
The spread in the values for the dynamic exponents obtained from the dispersion relations 
for small q and the median frequencies is considerable. In order to reconcile all estimates 
with Eq.( [4.1l|) almost the full width of the error interval (one standard deviation) is needed. 
However, Eq. fl4.lip represents a reasonable mean value of all estimates discussed so far and 
we therefore adopt it as our final estimate for z m and use z$ + z m = 3 to determine z^. 

As a final test of Eq.( |4.11| ) we use Eq. (fO|) in order to obtain a scaling plot of the dynamic 
structure factor. For qL = 2n, L = 30, 40 and 60 and z m = 1.38 the resulting scaling plot for 
S zz (q,uj) is shown in Figjl^. The corresponding plot for S xx (q, ui) with z^ = 3 — z m = 1.62 
is displayed in Fig.[13|. Data collapse can only be expected for sufficiently small us. Scaling 
works quite well for all frequencies up to the spin wave peak, where the intensity decreases 
slightly with increasing L, but the error bars are still overlapping. The quality of the data 
collapse in not so evident from S xx (q, lj) (see Fig.|l3|). At uj = the scaled intensity for 
L = 60 deviates from the ones for L = 30 and 40 by about one standard deviation. This 
could be an effect of the finite statistical sample. The line shape of the central peak in 
S xx still fits rather well to a Lorentzian. Only in the vicinity of u = are the data also 
compatible with the Gaussian lineshape expected from mode coupling theory [^2 . 



The above scaling analysis provides quite strong evidence for a violation of dynamic 
scaling in the sense that two different dynamic exponents are required in order to obtain 
scaling in the dispersion relations, the median frequencies, and the two components of the 
dynamic structure factor. However, we cannot prove from our data whether the estimate 
u w = zj, — z m = 0.24, which indicates the violation of dynamic scaling, corresponds to the 



transient exponent u w |2j,^J or if the measured difference consitutes an effective exponent 
for system sizes which are too small to ignore additional corrections. As pointed out in Sec. I 
the presence of energy conservation in our spin dynamics simulation gives rise to such cor- 
rections governed by the exponent ajv p3 |. On the basis of our data we cannot exclude the 



possibility, that these corrections yield the dominant contribution to the measured effective 
exponent u w . For a clarification of this point more information is needed from analytical 
theories of dynamic finite-size scaling [B3 and from simulations |BD . 



V. SPIN TRANSPORT AND THERMAL CONDUCTIVITY 

We have already mentioned in the introduction that the transport properties of the XY 
model near the critical point provide lattice analogues of the corresponding transport coef- 
ficients of 4 He near the A-transition. The conserved out-of-plane component M z of the mag- 
netization is the lattice analogue of the entropy density in 4 He and its associated transport 
coefficient therefore corresponds to the thermal conductivity of 4 He which is of experimental 
interest. In principle the thermal conductivity can be extracted from an extrapolation of 
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the characteristic frequency luq of S zz given by ujq 1 = \S zz {q,ti) |5|,|40| to q = . How- 
ever, in order to obtain a reliable extrapolation at the critical point a very high momentum 
resolution for q — ► is required and this can only be realized with very large systems (see 
Fig.|9] and Ref. ||40|| ). We therefore resort to an alternative approach already considered in 
Ref. P0| , we express the thermal conductivity by a current-current correlation function for 



a suitably chosen current density = (ji,k, j2,k, js,k) at lattice site k ||41|| . 

In order to identify the current density we reexamine the z-component of the equation 
of motion (see Eq. (|2.2| )) which reads 

isi = ~J E (SfSt-S?S* k ), (5.1) 
ai ieNN(k) 

where the sum is over all nearest neighbors of lattice site k. We define the ith component 



ji : k (i = 1, 2, 3) of the current density associated with the lattice point k by |40[ 

ji,k = J {^k^k+ei ~ ^k^k+e^) > (5-2) 

where the notation k + denotes the nearest neighbor of the lattice site k in the ith lattice 
direction. For the case of the simple cubic lattice studied here ei, e2, and can be visualized 
as the unit vectors of a cartesian coordinate system aligned with the lattice. The lattice 
divergence of the current density according to Eq.([5.2|) at lattice site k is then given by 



V-j* = = J E (SfSl-SfSI) (5.3) 

«=1 leNN(k) 

which is just the negative r.h.s. of Eq. (|5.1[). Note that the lattice spacing has been set to 
unity. In the spirit of hydrodynamics Eq. ( |5.1|) can be interpreted as second order discretiza- 
tion of an equation of continuity of the form dm^jdt = — V • j^. The density = S% and 
the current defined by Eq.( |5.2|) are located on staggered meshes, i.e., the current density 
component associated with lattice site k is located half way between k and k + ej, so 
that each of the finite differences appearing in Eq.( |5.3j ) is located at lattice site k just as 
the l.h.s. of Eq. (|5.1|) . This displacement of spins and currents on the lattice has no further 
consequences in our case, because we have applied periodic boundary conditions in all lattice 
directions. According to Eq.(44) of Ref. |^T[ we find the expression 



1 2 f°° 

D(q,u) = — — / V (ji,o(0)ji,i(t)) cosq- Ri cosut dt (5.4) 

k B TXzz 7T Jo Y 

for the transport coefficient D(q,u), where q = (q, 0, 0) and the same normalization as in 
Eq.( [4.2|) has been used. Note that the out-of-plane static susceptibility Xzz = (M^) / {ksTL 3 ) 
needed for normalization in Eq.( |5.4| ) is essentially constant as a function of system size L 
at fixed temperature. In the following we will analyze D(q,u) only for T = T c , where A = 
D(0,0) corresponds to the thermal conductivity measured in 4 He experiments. According 
to Eq. (|5.1| ) the current density has the scaling dimension 1 — z m — d/2, because Xzz has 
the scaling dimension zero. From Eq. ( |5.4|) we then find 2(1 — z m — d/2) + d + z m = 2 — z m as 
the scaling dimension of D(q, uj) so that naive finite-size scaling yields D(0, 0) = A ~ i, 2 -*™ 



at the critical point |25]. We therefore expect the scaling form 
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D(q, uu)k B T cXzz = L 2 - z -A(qL, uL z ™) (5.5) 

at T = T c , where ksTcXzz is just a normalization factor. A corresponding scaling plot of 
D(q,u)xzz versus u/JL Zm for L = 30, 40, and 60 and q = is shown in Fig.|i~4"|, where 
we have used our estimate z m = 1.38. The statistical noise in the data for D(q,u>) is 
considerably larger than the noise in the data for the structure factor (see also Ref. [|40|j ). 
The spread of the data points in Fig.|TJ| is of the same magnitude as the scatter of the data 
in each individual data set displayed. In view of these statistical uncertainties our data scale 
reasonably well for small uj and confirm the scaling exponent 2 — z m ~ 0.62 (see Eq.( [4.11| )) 
of the transport coefficient D(q,uj). For q ^ D(q,uj) scales accordingly as shown in Fig.[L5| 
for qL = 2n. The position of the spin wave maximum appears to be shifted to the right as 
compared to the spin wave peak in S zz (q,u>). According to Eq. ( |5.5|) scaling is obtained up 
to a value of about 0.6 of the scaling argument for both q = and qL = 2tt. According to 



Fig.[15] the spin wave maximum itself appears to be outside the scaling regime for the system 
sizes used here. For q ^ D(q,u) tends to zero as uj — > within the statistical errors. 

The thermal conductivity A = D(0, 0) is shown in Fig {II] as a function of the system size. 
The overall behavior of A with L is captured quite well by the expected power law (full line) 
for L > 30 and even the data point for L = 24 is only one standard deviation away. Systems 
with L = 20 or less may be too small to be in the scaling regime. In view of the considerable 
statistical error in the data slowly varying corrections to scaling as discussed in the previous 
section cannot be identified. In any case more theoretical information in needed in order 
to provide a reliable background for the interpretation of spin dynamics data of transport 



coefficients in the critical regime [38,39 . 



VI. SUMMARY AND CONCLUSIONS 

The easy-plane Heisenberg ferromagnet belongs to the XY universality class which has 
been demonstrated by the evaluation of the critical exponents and the static structure fac- 
tor. Unlike the standard plane rotator XY model the planar ferromagnet is endowed with 
a reversible spin dynamics which can be efficiently simulated by recently developed decom- 
position methods. Due to spatial and temporal symmetries the spin dynamics of the planar 
ferromagnet is expected to be in the same dynamic universality class as 4 He near the super- 
fluid normal transition, but may have different corrections to scaling. The data have been 
compared to field-theoretic and mode-coupling predictions with the following main results. 

1. Above the critical regime the spin dynamics data show a strong Lorentzian central 
peak in the in-plane component S xx (q,u>) in agreement with mode coupling theory. 
The out-of-plane component S zz (q, uj) displays a pseudo spin- wave peak which is also of 
Lorentzian shape and becomes increasingly prominent as q is increased. The presence 
of this peak is in accordance with field-theoretic predictions for T > T c . 

2. Below the critical regime strong spin wave peaks occur in both components of the 
dynamic structure factor. The shape of these peaks is captured very well by Lorentians 
as expected from mode-coupling theory. In addition to the spin wave peak S xx (q,uj) 
displays a central peak of Lorentzian shape and additional multi spin-wave peaks which 
do not appear in S zz (q, uj). For small q the dispersion relation u)(q) is linear in q and 
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the line widths T xx (q) and T zz (q) are quadratic in q as expected from mode-coupling 
theory. The qualitative agreement with the spin dynamics data for d = 2 suggests 
that the dynamics of the two-dimensional XY model may not be captured by vortex 
dynamics theories. 

3. At T c , in contrast to mode-coupling theory and in agreement with field-theoretic pre- 
dictions, S xx {q.uj) and S zz (q, uj) require different dynamic exponents in order to obtain 
scaling: z m = 1.38 ± 0.05 and = 3 — z m = 1.62 ± 0.05, whereas mode coupling 
theory yields z m = z^ = 3/2. The out-of-plane component S xx (q, uj) is dominated by 
a strong central peak. A shoulder at a finite frequency indicates the presence of a spin 
wave signal. In S zz (q, uj) a strong spin wave peak remains the dominating feature. 
The line shapes are still compatible with Lorentzians, the central peak in S xx is only 
compatible with a Gaussian shape very close to uj = 0. 

4. The transport coefficient D(q,uj), which provides access to the lattice analogue of 
the thermal conductivity of 4 He within the XY model, has been investigated at T c . 
The statistical fluctuations of the data are much stronger than those in the data for 
the structure factor which makes the scaling analysis less unique. However, scaling in 
agreement with the previously obtained dynamic exponents is obtained. The transport 
coefficient also shows a strong spin wave resonance for finite q and vanishes for uj — > 
in this case. The thermal conductivity A = D(0, 0) scales with the system size in the 
expected way within the error bars for sufficiently large systems. 
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FIGURES 

FIG. 1. Scaling plot for G;(q)/£ 2 ^ at T = T c for L = 20 (O), 40 (+), 60 (□), and 80 (x) as 
a function of qL with 77 = 0.035. Statistical errors of the data are much smaller than the symbol 
sizes. The dashed line displays the scaling function gi(qL) as given by Eq.( |3.8[ ) for a = 3.70 and 
b = 4.35. Note that qL = 20ir corresponds to the Brillouin zone boundary for L = 20. 

FIG. 2. In-plane component S xx (O) and out-of-plane component S zz (+) of the dynamic 
structure factor in the (100) direction for T = 1.1T C , L = 30, and q = ir/15 as functions of 
the dimensionless frequency uj/J. For small enough uj the line shapes are well approximated by 



Lorentzians given by Eq.(4.3) (solid line) and Eq.(Oh (dashed line), respectively. The ranges over 



which the Lorentzian approximation is valid differ significantly for S xx and S, 



zz ■ 



FIG. 3. Out-of-plane component S zz of the dynamic structure factor in the (100) direction for 
T = 1.1T C , L = 30, and q = vr/15 (O), q = 2vr/15 (+), q = 3vr/15 (□), q = 4vr/15 (x) as a function 
of the dimensionless frequency to/ J. 

FIG. 4. In-plane component S xx (O) and out-of-plane component S zz (+) of the dynamic 
structure factor in the (100) direction for T = 1.1T C , L = 30, and q = ir as functions of the 
dimensionless frequency u/J. A shoulder appears at the position of the spin wave peak in S zz . 

FIG. 5. In-plane component S xx (O) and out-of-plane component S zz (+) of the dynamic 
structure factor in the (100) direction for T = 0.9T C , L = 30, and q = ir/15 as functions of the 
dimensionless frequency u/J. Apart from the dominant spin wave peak additional resonances 
appear in S xx (arrows) at low intensities (see main text). 

FIG. 6. Dispersion relation uj(q) (O) in the (100) direction for T = 0.9T C and L = 30. The 
solid line is the T = dispersion relation obtained from linear spin wave theory (see Eq.( [4,6| )). 
The dashed line is a fit to the first two terms of the Fourier series given by Eq. Q4.7| ). 

FIG. 7. Line widths (a) T xx {q) and (b) T zz {q) of the spin wave peak in S xx for the (100) 
direction, T = 0.9T C , and L = 30. The solid lines are parabolic fits to the first four data points (see 
FigM). The dashed lines are fits with the first two terms of the Fourier series given by Eq.( 



FIG. 8. Extrapolation of the line widths T xx (q) (solid line) and T zz (q) (dashed line) to q = 0. 
The line widths vary as q 2 for small q as predicted by mode coupling theory (see main text). 



FIG. 9. Dispersion relation uJ zz (q) in the (100) direction for T = T c and L = 24 (+) and 40 
(O). The solid line displays a linear extrapolation of the data for q < ir/2 to q = 0. Deviations 
from linearity become visible only for q < it/ '20 where critical effects set in. 
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FIG. 10. Small q behavior of the dispersion relations (a) uj zz (q) and (b) u) xx (q) in the (100) 
direction for T = T c . The solid line is a power law fit to the data which yields z m = 1.38 ± 0.05. 
The dashed line is a power law with the exponent z^ = 3 — z m = 1.62 ± 0.05. The data for the 
smallest systems L = 20 and L = 24 deviate systematically from the power law. 

FIG. 11. Median frequencies u) xx (O, solid line) and cu^f (+, dashed line) for qL = 2ir. The 
solid and dashed lines display fits to power laws (see Eq.(|l|)). 



FIG. 12. Scaling plot of S zz (q, u)/ ((L/10) z ™G zz (q)) versus u/J(L/10) z ™ for qL = 2tt and 
system sizes L = 30, 40, and 60. 

FIG. 13. Scaling plot of S xx {q,u) / {{L/lQ) z <t>G xx {q)) versus U)/J{L/1Q) Z * for qL = 2vr and 
system sizes L = 30, 40, and 60. 



FIG. 14. Scaling plot of D(q, to)x zz /(L/W) 2 ~ Zm versus lu/J{L/W) z ™ for q = 0, L = 30, 40, 
and 60, and z m = 1.38. The value D(q = 0,lu = 0) corresponds to the thermal conductivity. 

FIG. 15. Scaling plot of D(q, u)x zz /(L/W) 2 - Zm versus uj/J(L/10) z ™ for qL = 2vr, L = 30, 40, 
and 60, and z m = 1.38. D{q,uj) shows a strong spin wave resonance and it vanishes for uj — > 0. 

FIG. 16. Thermal conductivity Xxzz versus system size L for L = 20, 24, 30, 40, and 60. The 
solid line displays the power law L 2 ~ Zm for z m = 1.38 for comparison. The power law represents 
the data reasonably well for L > 30. 
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